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ABSTRACT 

We present 230 realizations of a numerical model of planet formation in systems without gas giants. 
These represent a scenario in which protoplanets grow in a region of a circumstellar disk where water 
ice condenses and the surface density of solids is enhanced (the "ice line" ) , but fail to accrete massive 
gas envelopes before the gaseous disk is dispersed. Each simulation consists of a small number of 
gravitationally interacting oligarchs (protoplanets) and a much larger number of small bodies that 
represent the natal disk of planetesimals. Time zero of each simulation represents the epoch at which 
the gas has disappeared, and the dynamics are integrated for 5 billion years (Gyr). We investigate 
systems with varying initial number of oligarchs, oligarch spacing, location of the ice line, total mass 
in the ice line, and oligarch mean density. Systems become chaotic in ~ 1 Myr but settle into stable 
configurations in 10-100 Myr. We find: (1) runs consistently produce a 5-9 Af® planet at a semimajor 
axis of 0.25-0.6 times the position of the ice line, (2) the distribution of planets' orbital eccentricities 
is distinct from, and skewed toward lower values than the observed distribution of (giant) exoplanet 
orbits, (3) inner systems of two dominant planets (e.g.. Earth and Venus) are not stable or do not 
form because of the gravitational influence of the innermost icy planet. The planets predicted by our 
model are unlikely to be detected by current Doppler observations. Microlensing is currently sensitive 
to the most massive planets found in our simulations, and may have already found several analogs. 
A scenario where up to 60% of stars host systems such as those we simulate is consistent with all 
the available data. We predict that, if this scenario holds, the NASA Kepler spacecraft will detect 
about 120 planets by two or more transits over the course of its 3.5 yr mission. Furthermore, we 
predict detectable transit timing variations exceeding 20 min due to the presence of additional outer 
planets. Future microlensing surveys will detect ^ 130 analogs over a 5 yr survey, including a handful 
of multiple-planet systems. Finally, the Space Interferometry Mission (SIM-Lite) should be capable 
of detecting 96% of the innermost icy planets over the course of a 5 yr mission. 

Subject headings: celestial mechanics — planets and satellites: dynamical evolution and stability — 
planets and satellites: formation — planet-disk interactions — planetary systems 



1. INTRODUCTION 

The overwhelming majority of more than 400 exoplan- 
ets detected thus far are g as giants with masses compa- 
rable to Saturn or Jupiter (jCumming et al.|[2008D . Most 
methods of planet detection, including the Doppler radial 
velocity technique which has detected the overwhelm- 
ing majority of the known exoplanets, are biased to- 
ward higher planet mass as well as shorter orbital pe- 
riod ()Nelson fc Angellll998l: lUdrv et al.|[20Ql iCummind 
120041 ). Doppler surveys are seriously incomplete for 
semimajor axes larger than 5 AU, but an extrapola- 
tion of a debiased Doppler sample with a "flat" distribu- 
tion predicts that only ~17% of planetary systems con- 
tain g iant planets within 20 AU (jLineweaver fc Gretheii 
l2003t iCumming et al.l I2008D . By analyzing the sample 
of microlensing plane t detections i n a s urvey of high- 
magnification events, iGould et al.l (|2010( ) estimate that 
36% ± 15% of the host stars (with typical mass of '-^ 
0.5 Mq) host giant planets, with 0.02 Mj <M <h Mj 
per logarithmic decade in separation and mass. By com- 
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bining the results of ICumming et aLl (120081) fo r planets 
with a < 2 AU with those of lGould et al.l (|2010l ) for 2 AU 
< a < 20 AU, we estimate that 37% ± 13% of stars host 
giant planets (m > 0.3 Mj). We conclude that the frac- 
tion of stars hosting giant planets with a < 20 AU is 
likely to be at least 20% but less than 50%. 

Thus a substantial fraction, and probably the major- 
ity of stars do not host giant planets within < 20 AU. 
However, studies of star-forming regions of different ages 
have shown that all or nearly all solar-mass stars be- 
th eir lives with disks (iHaisch et al.ll2001l: iLada et al.l 
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20061: iSicilia-Aguilar et al.l l2006t iLuhman et al.l 120081 
20101 : iMassi et al.ll2010f l. If the accretion of solids is suf- 
ficiently rapid and efficient in these disks, then plane- 
tary systems are presumably equally numerous around 
middle-aged stars and the relative paucity of gas giants 
demands explication. This conclusion is also supported 
by th e limited statistics of debris disks aro und solar mass 
stars (|Wvattll2008l : iCarpenter et al.l[2009l ). 

A simple explanation is that gas giants never form 
around the majority of stars. The core accretion theory 
of giant planet formation predicts this outcome if disk 
gas usually disperses before the growth of a sufficiently 
massive solid core triggers runaway accretion of the gas. 
The th reshold mass is curre ntly thought to be at least 
5 Mq (jHubickyj et al.l[2005| ) but depends sensitively on 
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the gas opacity (jAvliffe fc Batel[2009[ ). Canonical models 
of accretion in a minimum-mass Solar Nebula (MMSN) 
fail to produce a sufficiently massive core at the orbit of 
Jupiter (Polla ck et al.l 119961 ) in the ~2-6 Myr timescale 
on which disks are o bserved to dissipate (H aisch et alJ 
120011: iHubickvj et alJ[2bO& lEvans ct al. 20091 . A com- 
mon explanation for the formation of Jupiter is to 
include an "ice line" at 3-5 AU where water con- 
denses and the surface density of s olids is substantially 
elevated, promoting core growth (Stevenson & Lunind 
1988; P ollack et al. 1996; Kokubo & Ida 2002; Ida & Lin 
2004bt ). Disks around more massive and/or more metal- 
rich stars presumably hav e greater amounts of solids, 
accelerating core accretion ([Laughlin et al.l [2004 iCurrid 
I2009D , consistent with the observed correlation of giant 
planet frequenc y with host star mass a nd metallicity 
jGonzalcz 1999; 'Fischer fc Valentil l2005t I Johnson et"all 
[2007; Kcn nedv fc Kenvon 2008b[ ). Theoretical studies 
also predict a correlation between disk surface density 
and the frequenc y (and rnass) of giant planets and lower 
mass icy planets (|Rafikovl200l:lKenvon fc Bronilev[|2009l 
I20T0I) . 

A second explanation is that giant planet formation 
may be stymied by the tendency of cores to migrate 
inwards wh ere there is insuf ficient gas to form a gi- 
ant planet (jlda fc LinI l2008bD . This scenario is pred- 
icated on the operation of Type I migration in which 
torques from the gas disk become important for bodies 
more massive than Mars. The planet-metallicity correla- 
tion is explained if cores in disks with more solids grow 
more rapidly to the threshold of runaway gas accretion, 
thereby opening up a gap in the disk and halting Type 
I (but not Type II) migration. However, the magnitude 
and eve n the si gn of Type I migration remains v ery un- 
certa in (|Li et a l. 2009; Muto fc Inutsuka 2009; Y u et al.l 

Another possibility is that giant planets are ubiq- 
uitous but have migrated or been s c attered outward 
(ICrida et al.l [200I IVeras et al.l [2001 IScharf fc Mel^ 
120091 ) to distances where they are detectable only 
by microlcnsing when they b e have like isolate d 
lenses (Pi Stefano fc Scalzdll999"al[bl : IGould et al.ll2nn7n . 
or by their infrared emission under exceptional 

circums tances of youth, r nass, and proxirn ity to 

Earth (iDebes fc Si gurdsso ij [2001 iKalas et al.l [20081: 
iMarois et al.l 120081 : iThalmann et al.l l2009( ). Finallv. gi- 
ant planets may m igrate inwards to disrup tion within 
the Roche zone (jPatzold fc Raueil 120021 ) . although 
there are limits on the ubiquity of s uch occurrences 
(jPinsonneault et al.|[200lHQuillen|[200l . 

The first two explanations predict that systems lack- 
ing gas giants will contain "failed" cores of Earth 
to Neptune mass tha. t prefe r entially formed near 
the ice line (llda fc LinI l2004al: iThommes et all 120081 : 



iKennedv fc Kenvonll2008bt iMordasini et al.ll2009l ). Un- 
less such objects migrated inward, they would remain 
invisible to the Doppler technique: the signal from a 
10 M0 planet at 5 AU is 0.6 m s~^, well below the sta- 
bility of radial velocity mea surements on decadal base- 
lines (jCumming et al.ll2008| ). However, the gravitational 
microlensing technique is capable of uncovering such 
planets and sever al have already been found. Indeed, 
iSumi et"aL[ (|2010[ ) argue that the slope of the mass (ra- 
tio) function for planets beyond the ice line is quite steep. 



such that Neptune-mass planets are ~ 7 times more 
common than Jup iter-mass planets. Combined with the 
IGould et al.l (|2010[ ) normalization of the frequency of gas 
giants in this region, this implies that the majority of 
stars host Neptune or lower mass planets beyond the 
snow line. 

In this paper we investigate a scenario in which sev- 
eral protoplanets or "oligarchs," but not giant planets, 
have formed at or beyond the ice line at the time disk 
gas has disappeared. A major premise of our initial 
conditions is that Type I migration was not effective in 
these disks. (Type II migration does not act on planets 
much less massive than Jupiter) . We carry out direct nu- 
merical integrations of the orbital and mass evolution of 
these protoplanets as they accrete additional mass from 
a residual disk of planetesimals. We integrate the orbits 
of the (proto)planets over 5 Gyr to ascertain the stabil- 
ity of these systems and their configuration at a plausible 
epoch at which they might be observed. 

Our si mulations c omplement the works of 'Ida fc Lira 
(|2008al) : IKennedv fc Kenvon (2008a); Mordasin i et al.l 
(^2009"), and 'Thommcs ct al. (2008). Ida fc LinI (|2008aD 
and Mordasini ct al. ( 2009) use analytical models of the 
orderly growth of cores in disks during the interval that 
gas is present. They include migration due to torques ex- 
erted by the disk, but neglect subsequent, chaotic grav- 
itational interactions between the core s and the resid- 
ual di sk of solids (Ida fc Lin 2008b). IThommes et al.l 
(12008) model the dynamical interactions between pro- 
toplanets/cores as well as the gas disk but analytically 
proscribe accretion from a fixed disk of solids. In cases 
of low disk mass or rapid gas removal, all three investiga- 
tions pred ict formation of Earth- to N eptune-sized bod- 
ies. Like IKennedv fc Kenvoiil (|2008a[ ). our simulations 
begin at the end of the orderly growth phase when the 
gas has disappeared, and assume that no giant planets 
have formed, but unlike them, we assume that a massive 
residual disk of smaller bodies is still present. 

Our work also contrasts with investigations of the 
dynamical evolution and configuration of systems of 
giant planets like t hose d e tected by Doppler surveys, 
e.g.. iRasio fc FordI (Il996(). lAdams fc Laughlii] ([2005) . 



Chatteriee et all (120081 ) . lllavmond et al.l ()2009a[) . and 
Ravmond et al.l (|2010[ ). We expect that the evolution 






of systems of solid planets that emerge from the icy part 
of a disk will differ substantially, primarily because the 
Safronov number 

S^^ (1) 

will be < 1, whereas S* ^ 1 in systems of giant planets 
(excluding "hot" Jupiters) . More efficient accretion, less 
intense scattering, and stronger coupling to the planetes- 
imal disk are expected. 

The goal of our simulations is threefold: first, we want 
to predict the evolution and final configuration of such 
systems for a range of plausible initial conditions. Sec- 
ond, we wish to determine if and how such planets could 
be detected by present or future means, i.e., the Ke- 
pler and Space Interferometry Missions (SIM-Lite) and 
ground-based microlensing surveys. Lastly, we want to 
establish how measurements of these objects might be 
used to infer the initial conditions and histories of these 
systems, especially important given how poorly planet 
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Fig. 1. — Schematic of the typical initial configuration of our 
systems. Oligarchs are shown in gray in a zone of enhanced surface 
density between 5 and 6 AU. Small bodies are shown in black 
and represent a much larger number of planetesimals in the disk. 
Plotted along the bottom of the graph is an approximation of the 
logarithm of the surface density of material. Beyond the region of 
oligarch formation, the density of planetesimals follows a S <x a~^ 
mass distribution. The inner system is left empty because this 
region has only a minor effect on the evolution of the outer system, 
but requires considerably more time to simulate. 

formation is understood. 

2. METHODS AND MODELS 
2.1. Approach and Assumptions 

Each numerical realization consists of a IMq cen- 
tral star surrounded by a planet-forming disk (Fig- 
ure [1]). Our scenario assumes the canonical theory of 
planet formation which consists of three phases: (1) run- 
away accretion of protoplanets from a disk of planetes- 
imals; (2) slower oligarchic growth of these protoplan- 
ets as they consume neighboring planetesimals and each 
other; and (3) a chaotic or giant impact phase when 
the mass in residual planetesimals falls below that in 
the protoplanets and the oligarchs' orbits begin to cross 
(JGoldreich et al. 2004; Kcnyon & Bromley 2006). The 
disk includes a region of width 6ice immediately beyond 
the ice line (oice) in which the surface density of solids is 
enhanced by the transport of water vapo r out of the inner 
disk and condensation as ice at a > ajce (jCuzzi fc Zahnlg 
[2004; Ciesla & Cuzzi 2006). We assume that oligarchic 
protoplanets have appeared in this region because of the 
enhanced density and the relatively short orbital time 
scale compared to the disk further out. A background 
disk of unincorporated planetesimals extends from Ojce 
to 400 AU following the surface de nsity profile S oc a~^ 
found in observations of protostars ([Andrews fc Williamsl 
I2007f ). For computational efficiency, we exclude the re- 
gion outside of 100 AU and the region inside the ice line. 
A few sets of simulations included inner bodies to as- 
say their effect on the formation of planets further out 
(Section [m. 

We assume a disk of solar composition with the total 



Fig. 2. — Fraction of mass in small bodies within the orbit of 
the outermost oligarch during the first 1 Gyr for a run from the 
ST3 set. More than 95% of the mass is in the small bodies at the 
start of the simulation. By 500 Myr the region has been almost 
completely cleared of small bodies. Only as the outer oligarch 
migrates outward do more small bodies enter the region. 

mass within 400 AU of 0.03 Mq and a corresponding 
mass of conde nsible solids (rock and ice) of ~ 150 M^ 
([Loddersll2003n . The mass of the background disk within 
the 100 AU simulation region is 43-48 M(^, depending 
on how much mass is moved into the ice line. The addi- 
tional mass added to the ice line is varied (Section 12. 2p . 
We specify the number of oligarchs n and the spacing 
between them in Hill radii 6; this sets the mass of each 
oligarch and the total mass in oligarchs. The remaining 
mass, in fact the majority of the mass in all simulations, 
is distributed evenly among small bodies that represent 
primordial planetesimals. The number of small bodies is 
limited by computational resources and is 500 in all of 
our simulations, except for a single run with 1000. The 
mass of each small body is about O.OSM^, much less than 
the oligarch masses. We do not model the fragmentation 
of planetesimals, the production of dust by a collisional 
casc ade, and the removal of that dust by stellar radia- 
tion ([Wvattll2008f ). We discuss the possible consequences 
of fragmentation on our conclusions in Section 15.31 Oli- 
garchs and small bodies were given non-zero random in- 
clinations (|z| < 10"'^ and 10~^ degrees, respectively) 
and eccentricities (e < 10~^ and 10~^, respectively), al- 
though we experiment with higher initial i and e values 
in a single simulation set (see Section [^75]) . 

The equations of motion of each particle were inte- 
grated using the hybrid integrator code Mercury6 with 
the combination of a second-order, mixed-variable sym- 
plectic integrator and the Burlirsch-Stoer integrator for 
close encounters (|Chanibers 1999). Mercury6 divides the 
simulation into large and small bodies. Large bodies in- 
teract gravitationally with, and can collide with both 
large and small bodies. Small bodies also interact gravi- 
tationally and collide with large bodies, but they cannot 
collide with other small bodies. Computations were per- 
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formed on the TeraGrid network (jCatlett et al.ll2007[ ). 
Integrations were performed for 5 Gyr, except for more 
computationally intensive runs used to check our condi- 
tions (Section 12. 3|) . Ten replicate simulations were per- 
formed for each set of parameters. The initial positions 
and orbits were varied slightly between simulations in 
a set. A timestep of 40 days was used for most sim- 
ulations, based on the expected orbital period of the 
innermost oligarch and the requirement that there be 
at least 10-20 time steps per orbital period, a conser- 
vative setting (Raymond ct al. 2010). Simulations that 
included oligarchs in the inner system had a timestep of 
8 days, and simulations with a close-in ice line were given 
a timestep of 20 days. Runs where there were fewer than 
10 timesteps per oligarch orbit for an extended period of 
time were either adjusted, rerun, or had short test sim- 
ulations run parallel to them to test the accuracy of the 
results. 

For computational efficiency, small bodies in most runs 
are removed after 1 Gyr. Previous work has shown 
that this economy will not significantly affect the evo- 
lution of the oligarchs if the mass surface density of 
the small bodies is much less than that of the oligarchs 
(jKenyon fc BromlevI 120061 ). Figure [2] shows the fraction 
of mass in small bodies inside the orbit of the outermost 
oligarch for the first Gyr of a run in our standard star. 
Although at the start of the simulation the small bodies 
represent the bulk of the mass, by 100 Myr they are less 
than 50% of the mass in the oligarch region. By 1 Gyr, 
small bodies account for only ~ 10% of the total mass. 
We do not observe significant changes in the orbits of the 
oligarchs as a result of the removal of the small bodies 
at 1 Gyr. Oligarchs tend to move into stable orbits long 
before 1 Gyr. This can be seen in Figures [3] and H) In a 
single run from each set we retain the small bodies for 2 
Gyr as a check. 

2.2. Parameter Values and Initial Conditions 

Our simulations are described by four principal param- 
eters: Cice, -Mice I ^i ^nd n. We also vary the mass density 
of the oligarchs p. Table [T] lists the parameters for our 
18 sets, each of which consists of 10 replicate runs. Five 
additional sets of 10 sensitivity simulations are shown at 
the bottom of the table and discussed in Section 12.31 



Ice line location: The location of the ice line Ojce de- 
pends on the opacity and rate of viscous dissipation in 
the disk and the mass of the central star, and is time 
dependent (jCiesla fc Cuzzil |2006|) . The ice line in the 
primordial solar system has been variously placed near 
5 AU (to stimulate the rapid formation of Jupiter's core 
and explain its icy satellites) ([Stevenson fc Luninelll988| ) 
or at 2.7 AU (to coincide with the transition betwe en hy- 
drated and anhydrous asteroids) ()Abe et al.l 1200(1). The 
position of the ice line presumably varies between plan- 
etary disks. In the majority of our simulations, we set 
Oice — 5 AU. In the CI sets (30 runs) a^^e was fixed at 
2.7 AU. 

Ice line mass: Our choice of mass just beyond the ice 
line is guided by an estimate of the mass of water trans- 
ported through the inner disk that re-condensed at the 
ice line, and the mass of solids in the cores of the outer 
planets in our solar system. To the small amount (< 1 
M^) of solids within Oice < a < ai^e + ^ice predicted by 
a simple E oc a~^ distribution, we add a fraction of the 
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Fig. 3. — Evolution of the system in a run from the ST4 set. 
Oligarch swaps (exchange of order with distance from the star) 
are marked with dashed lines, excluding the first 10 Myr when 
oligarch swaps are more frequent. Because oligarchs often undergo 
numerous swaps when they approach each other's Hill radius, each 
dashed line may represent multiple swaps. 



0.4 










E 


0.3 

0) o.a 










; 


0.1 




\!l>k I ' 


l\:\i.k-:^ 


,-i;,;N tl ;tJf)l. ,1, 




0.0 






. , ' 


' V ■ '■,, " 1,, 


• ' . 1 11 _l , ; 


12 


- 








- 


10 


'- 








-_ 


S a 

< 

^^ 6 


f 


iJVw"*- 


y-'^V. 


— -_^ 


^ 


to 

4 


n^ 


^^ 






: 


2 


'r 








: 











: 












- 


6 


■^ 








- 






- 


S! 

a 

^ 4 


' " 




ID 

6 z 






^ 





- 








- 



200 



400 600 

time (Myr) 



BOO 



1000 



Fig. 4. — Evolution of the system in a run from the ST3 set 
showing the migration of the innermost migratory planet (IMP) 
(green). Initially, the IMP (green) migrate inward, while the sec- 
ond and third oligarch (blue and teal respectively) migrate outward 
through their exchange of angular momentum. All three oligarchs 
are driven inward at times by their interactions with the planetes- 
imal disk. After ~ 100 Myr the region from to 6 AU has been 
almost completely cleared of small bodies, causing the inner two 
oligarchs to settle into relatively stable orbits. This type of angular 
momentum exchange is common in simulations that start with 3-4 
oligarchs. 
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TABLE 1 

Simulation Initial Conditions 



Name 


6 


n 


A^ice 


a:™ 


rui 


^Comments 




[Rh] 


1 


(Me)(AU) 


(M®) 




ST2 


8 


2 


35 


5.0 


0.44 


Standard 


ST3 


8 


3 


35 


5.0 


0.44 


Standard 


ST4 


8 


4 


35 


5.0 


0.44 


Standard 


LM2 


8 


2 


10 


5.0 


0.44 


Low ice line mass 


LM3 


8 


3 


10 


5.0 


0.44 


Low ice line mass 


LM4 


8 


4 


10 


5.0 


0.44 


Low ice line mass 


CI2 


8 


2 


35 


2.7 


0.44 


Close ice line 


CIS 


8 


3 


35 


2.7 


0.44 


Close ice line 


CI4 


8 


4 


35 


2.7 


0.44 


Close ice line 


MP2 


5 


2 


35 


5.0 


1.78 


"Medium packed" oligarchy 


MPS 


5 


3 


35 


5.0 


1.78 


"Medium packed" oligarchy 


MP4 


5 


4 


35 


5.0 


1.78 


"Medium packed" oligarchy 


HP7 


2 


7 


35 


5.0 


0.70 


"Highly packed" oligarchy 


HP9 


2 


9 


35 


5.0 


1.98 


"Highly packed" oligarchy 


HP12 


2 


12 


35 


5.0 


3.63 


"Highly packed" oligarchy 


LD2 


8 


2 


35 


5.0 


0.44 


Low density oligarchs 


LD3 


8 


3 


35 


5.0 


0.44 


Low density oligarchs 


LD4 


8 


4 


35 


5.0 


0.44 


Low density oligarchs 



RD3'' 


~ 8 


3 


10 


5 


0.30-0.70 


TS3 


8 


3 


35 


5 


0.44 


EV3'' 


8 


3 


35 


5 


0.44 


103" 


8 


3 


35 


5 


0.44 


HE3'^ 


8 


3 


35 


5 


0.44 



Random m, a of oligarchs 
1000 small bodies 
Earth, Venus included 
Oligarchs in inner system 
Higher initial e, i 

[^fThe mass of the oligarchs in this system varies randomly by 10% and 
the semimajor axis by 0.2 AU . b was allowed to vary based on the 
mass and semimajor axis, but it was not allowed to go below b — 6 or 
above 6—9. 

^frhis simulation set was run to 100 Myr with a planetesimal disk and 
then to 1 Gyr without small bodies. 

[jThis system has 130 oligarchs in the inner system following a fe — 8 
and S oc a~ distribution. The total mass of the inner system was 
2.2 Mq) spread between 0.5 and 5 AU . This simulation set was run to 
100 Myr. 

[jThis system has small bodies with \i\ < 10^ and e < 5 X 10^ for 
the small bodies in the region of the oligarchs. This simulation set 
was run to 250 Myr. 



total amount of water from the disk outside this region. 
We adopt a minimum value of M^^^ = 10 Mq. Our max- 
imum value (Mice = 35 M^) assumes all four outer solar 
system planets formed near the ice line (jThommes et alj 
1200 2') and sums their initi al core masses: 10 M m for 
Jupiter, 15 M® for Saturn (|Hubbard et al.ll2009f ) and 5 
Mq (below the critical threshold) each for Uranus and 
Neptune. Mice — 10 Af© and 35 M© correspond to 
15% and 60% of the disk's water in the ice line and a 
background disk mass within 100 AU of 43 and 48 M©, 
respectively. Other estimates of the mass in the ice 
line are similar (IStevenson fc Luninelll988HKornet et alJ 
l200llIdafcLinll2008bn . 

Oligarch spacing: The spacing between oligarchs is 
specified as a multiple of the Hill radius, 



Rh = a 



3AL 



1/3 



(2) 



where m is the in itial oligar c h ma ss, and M* — 1 Mq. 
We use 6 = 8 fChambcrs' '2006') as well as fo = 5 
(|Ford fc Chiang 2007; Ravmond ct al. 2009b, 2010). Ad- 
ditional runs contain "overpacked" oligarchies with b — 2 
(the approximate Roche limit). This last value is well 
within th e instability lim it h = 3^' ^ where accretion can 
be rapid ( jGlad man"1993') but outside horseshoe (1:1 res- 
onance) orbits (|Collins fc Sari„2009, ). 



Numbers and masses of oligarchs: We vary the num- 
ber of oligarchs n between two and four in the 6 = 5 
and & = 8 simulation sets, in analogy to the number of 
cores that formed at the ice line in our solar system. For 
the case of n = 3 and a^ce = 5 AU, we fixed the width 
of the ice line 6\rp to 1 AU (jStevenson fc Lunind 119881 : 
iKornet et al.ll20(M ). The initial mass of each oligarch is 
related to n. 



^ice, and Sice by 



3M, f % 
\nba 



(3) 



This gives initial oligarch masses of 0.44 and 1.78 M© 
for 6 = 8 and 6 = 5, respectively. We subsequently scale 
(5ice with ajce and n so that for a given value of 6, the 
initial oligarch mass is unchanged. In the overpacked 
(6 = 2) scenario, the requirement that the core mass not 
exceed 10 M© (and seed giant planet formation) requires 
n > 4. In these case, we fix the total mass in oligarchs 
nm as 0.25, 0.50, and 0.75 times the total ice line mass. 
We calculate n using Equation ([3|). This gives n = 12, 
9 and 7 and oligarch masses of 0.7, 1.98, and 3.63 M©, 
respectively. In any given run, all initial oligarch masses 
are the same, except for the RD3 simulation set, in which 
the oligarch mass is allowed to vary slightly. 

Oligarch mass density : We u se mean densities p pre- 
dicted by IGrasset et al.l (|2009( ) for planets of 60% ice 
(15% for the LM sets) and the rest rock. However, super- 
Earth-mass bodies may retain significant envelopes of 
H/He gas, giving them mean densities more akin to that 
of Neptune {p = 1.64 g cm~^), and this would increase 
the cross-section for accretion. To investigate this effect, 
the LD simulation sets were run with p = 1 g cm^'^. 

2.3. Sensitivity Runs 

Larger number of small bodies: The 500 small bodies 
in each simulation represent a much larger population 
of planetesimals in the disk. One run (the TP3 set) 
includes 1000 small bodies, each with half the mass of 
those in the 500-small body runs. Neither 500 nor 1000 
small bodies is physical (there may actually be trillions 
of planetesimals). The goal here is to ascertain if the re- 
sults depended sensitively on the number of small bodies 
or their mass. 

Retaining small bodies for 2 Gyr: For a randomly se- 
lected run in each set (excluding the CI and HP sets), the 
small bodies are kept in the simulation for 2 Gyr instead 
of 1 Gyr. This is done to verify that removal of the small 
bodies does not create a bias in the results. 

Non-identical oligarch masses: A real system of proto- 
planets will not have identically spaced, identical mass 
oligarchs. In one set of sensitivity runs (RD3 set) we 
vary the mass by ^10% and the semimajor axis by ±0.2 
AU. In these simulations 6 is allowed to vary between 7 
and 9 as a result of the randomized semimajor axis. 

Mass in the inner system: Two sets of runs investigate 
the effect of mass in the inner system (a < ajce) on planet 
formation in the outer system. In one set (103), we in- 
clude a disk of 130 small (0.002-0.1 A/©) oligarchs with 
a mass surface density E ^ a^^, separation of 6 = 8, 
and a total mass of 2.2 M© (the total mass of the in- 
ner planets in the solar system). Given that a common 
outcome of higher-resolution A'^-body simulations of ter- 
restrial planet accretion are two planets of roughly equal 
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mass (jRavniond et alj |2009(J) . a second set of runs was 
performed that included two planets with the masses and 
orbits of Earth and Venus. 

Small Body Eccentricities: Our small bodies start with 
non-zero but very small values of orbital eccentricity and 
inclination. However, gravitational perturbations by the 
oligarchs will drive these values t o finite values (^ 0.05) 
even in the presence of disk gas (jChambersI I2006D . We 
run a single set of simulations to 250 Myr with |i| < 10^^ 
and e < 5 X 10~^ for the small bodies in the region of 
the oligarchs. 

3. RESULTS 

3.1. Dynamical Evolution 

We are primarily interested in describing the gross dy- 
namical evolution of these systems over the first 1 Gyr, 
especially when small bodies remain in the system and 
are being scattered or accreted by the oligarchs. At later 
times, most (but not all) of the systems do not substan- 
tially evolve. We analyze orbital parameters with a res- 
olution of 100 kyr: important events happen on shorter 
timescales than this, especially in the first 1 Myr of the 
simulations when systems are most chaotic. However, 
analysis of such events is not the focus of this work. 



TABLE 2 

Evolution Statistics 



Name 


Tclcar^ 


Resonance 


Frac Order 


Inner Planet 




(Myr) 


Crossings 


Preserved 


Migration (AU) 


ST2 


290 


100 


0.6 


3.7 


ST3 


80 


470 


0.2 


3.8 


ST4 


90 


980 


0.2 


3.7 


LM2 


450 


260 


0.7 


3.2 


LM3 


290 


820 


0.1 


3.0 


LM4 


240 


1490 


0.1 


3.2 


CI2 


62 


50 


0.4 


2.0 


CI3 


65 


320 


0.2 


1.9 


CI4 


56 


390 


0.3 


2.0 


MP2 


80 


70 


0.7 


3.8 


MPS 


43 


500 


0.6 


3.8 


MP4 


36 


1297 


0.2 


3.8 


HPT 


16 


3900 


0.1 


3.1 


HP9 


26 


11600 


0.0 


3.3 


HP12 


66 


28600 


0.0 


3.5 


LD2 


150 


80 


0.6 


3.8 


LD3 


74 


600 


0.3 


3.9 


LD4 


77 


1050 


0.3 


3.8 



ffi 



10 AU. 



i the average time for tlie disl^ to lose 70% of it's mass witliin 



We use five metrics to describe the evolution of each 
system (Table [5]): (1) the time Tdear in which 70% of 
small bodies are removed from inside 10 AU; (2) the 
number of mean motion resonance (MMR) crossings ex- 
perienced by the oligarchs with other oligarchs over the 
first 1 Gyr; (3) the fraction of runs in which two or more 
oligarchs "swap", i.e., exchange order with distance from 
the star; (4) the distance of inward migration by the (ul- 
timately) innermost protoplanet; and (5) the number of 
oligarchs ejected from a system. 

Disk clearing: r^iear is the time over which 70% of 
small bodies within 10 AU are accreted or ejected. The 
timescale is not sensitive to the precise choice of outer 
boundary. Tdear can be as long as several hundred Myr 
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Fig. 5. — Number of MMR crossings in 100 Myr bins, normalized 
by the first bin. This is a measure of the relative levels of chaos in 
the system. For our purposes, a resonance crossing occurs when- 
ever an oligarch crosses a 1:1, 2:1, 3:2, 3:1, 4:1, 5:3, or 5:2 mean 
motion commensurability with another oligarch. Total number of 
resonance crossings are reported in Table |51 



(Tabled . r^iear is shorter in systems with more oligarchs 
(more accreting bodies), lower values of b (oligarchs scat- 
ter each other onto more eccentric orbits); a closer ice 
line (shorter orbital period and dynamical time scale), 
and lower oligarch mass density (greater cross section of 
accretion). The LM sets have longer Tricar values be- 
cause there is less concentration of mass in the ice line, 
proximal to the oligarchs. The oligarchs do not grow as 
quickly nor scatter as efficiently. 

Mean motion resonance crossing: We count the num- 
ber of times an oligarch passes through a 1:1, 2:1, 3:2, 3:1, 
4:1, 5:3, or 5:2 mean-motion commensurability with an- 
other oligarch (Table [2]). The greatest number of MMR 
crossings occurs in runs with lower values of b and larger 
n. Figure [S] shows the normalized rate of MMR crossings 
per time for the six primary simulation groups, binned 
in 100 Myr intervals. In all cases the rate decreases with 
time and by 1-3 orders of magnitude over the first 1 Gyr 
as the systems evolve. 

Oligarch swapping: In 43% of the runs in the ST sets, 
at least one swap (where two oligarchs exchange rank 
in semimajor axis) occurs between 50 Myr and 1 Gyr. 
Swapping occurs when two oligarchs approach within two 
Hill radii (the zone of strong scattering). The oligarchs 
can collide, or they can ente r horseshoe orbits (within 
a single Hill radius) (jCollins fc Sarii.2009f ). In the latter 
case, they often exchange places quickly, usually in ^ 1 
Myr. Figure [3] shows an extreme case where there are at 
least 13 distinct swaps over the first 1 Gyr, excluding the 
chaotic period in the first 50 Myr. 

Inner Planet Migration: In all of our primary runs, an 
oligarch migrates inward to a position between 1 and 3 
AU, most often settling between 1.2 and 1.9 AU, i.e. 3- 
4 AU from its initial starting place. In 81% of cases this 
body has grown to become the most massive planet in 
the system. Most runs resemble that of Figure H] Fig- 
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ure [3] shows an unusual run in which two inner oligarchs 
migrate inward. The inward migration of one or more oli- 
garchs is a manifestation of the redistribution of angular 
momentu m in a circums tellar disk and its resulting radial 
dispersal (Pringle 1981). The angular momentum of the 
inwardly migrating oligarch(s) is lost to scattered small 
bodies and some of it is transferred to outwardly migrat- 
ing oligarchs (Figure |4]). This process will occur as long 
as there is a sufficiently massive disk of small bodies, i.e., 
for tens or hundreds of Myr (Table [2]). Analogous events 
may have unfolded during the early dynamical evolution 
of the outer Solar System: Jupiter migrated inward while 
the other giant planets moved outward as a result of an- 
gular moment um exchange through a residual disk of 
planetesim als CMal hotral Il993t iHahn fc Malhotral 119991: 
iGomes et a l. 2005). In our simulations, the Safronov 
number is less than or not much greater than one, and 
si gnificant accreti on of mass occurs during migration. 

Ilda et ail ()2000fl formulated the migration rate of a low 
mass planet moving through a planetesimal disk as 



da a AnTiptt^ 



dt 



K 



M, 



(4) 



where Pr- is the Keplerian orbital period and M» is the 
mass of the central star. The migration rate is indepen- 
dent of planet mass. Recast in terms of the ice line mass, 
the migration timescale is 



'migrate 



Pr 



2a M* 



(5) 



which as short as ~ 1 Myr for the undepleted disk. The 
observed migration timescale is slower (~10 Myr) and is 
probably in part due to the depletion of the disk by the 
oligarchs themselves, but may also reflect the inability of 
our simulations with low numbers of particles to correctly 
resolve the distribution of planetesimals in horseshoe or- 
bits that most strongly interact with the planet. 

Oligarch Ejection: Ejection of oligarchs can o ccur dur- 
ing the final, chaotic phase of planet formation ([Lissaueri 
[19871 IStevensonI [19991: IDebes fc Sieurdssonl [20071) . An 
oligarch is considered "ejected" in our simulations if it 
attains a > 400 AU. Ejection of an oligarch is frequent 
(13 of 30 runs) in the b — 2 simulations (Figure [5]). In 
most of these cases, an ejection occurs within 1 Myr af- 
ter two oligarchs appear to enter a resonance. Resonance 
between two oligarchs increases their orbital eccentrici- 
ties, and also excites neighboring oligarchs, and this is 
sometimes suffcient to eject smaller oligarchs from the 
system. In a single simulation, two oligarchs stayed near 
resonance for 100 Myr, causing the ejection of 3 other 
oligarchs. 

3.2. Configuration at 5 Gyr 

Mass: Figure [7| shows the final system configurations 
produced by 9 sets of 10 replicate runs (standard, MP, 
and HP sets) after 5 Gyr. Only 4 of the systems contain a 
single planet; all started with only n = 2 oligarchs. Mass 
segregation (the tendency of higher mass planets to ap- 
pear closer to the star) occurs under all conditions (Fig- 
ure[51). Of the 180 primary runs, ^ 94% produce a planet 
that ultimately resides between 0.25 and 0.6 a^cc, and in 
78% of all 6 = 8 runs, this planet is the most massive one. 



This effect is least pronounced in the 6 = 2 set. Approxi- 
mately one-fourth of the total mass in the ice line is incor- 
porated into these planets. In fewer than half of the runs 
did the initial innermost oligarch beco me this innermost 
planet (Figure [T]), in agreement with IChatteriee et al.l 
(f2008;i) , who found that in systems of equal- mass gas gi- 
ants, each planet has roughly equal probability of becom- 
ing the innermost. Planet mass decreases with semima- 
jor axis between 0.6 to 2.5 a^cc- Planets ending outside 
~ 2.5 fljce have masses close to that of the original oli- 
garch. These bodies were scattered outside the ice line 
early in the run and have accreted little mass. 




Fig. 6. — Evolution of the system in a run from the HP12 set 
showing the ejection of two oligarchs (teal and orange bodies) at 
~ 4.25 Gyr and ~ 4.5 Gyr. 



Eccentricity: We find no correlation between orbital 
eccentricity and mass, semimajor axis, or number of oli- 
garchs in our runs. Our runs with b — 2 (and the most 
oligarchs) produce a larger dispersion in eccentricity, but 
this could be a result of the higher total mass in oli- 
garchs. For & = 8, no planet ended with e > 0.3, 
and for the ST sets only 2 planets had eccentricities 
above 0.2. This contrasts with the observed distribution 
amongst detected planets, and simulations of systems of 
^s giants (Udrv & Santos 2007; Chattcriec ct al. 2003; 
IJuric fc Tremaind 12008) (Figure (H]). Smaller values of 
eccentricity are expected in systems where the mass of 
the oligarchs that perturb each other is lower relative to 
the mass of th e disk of small bodies that dampen such 
perturbations ([Chamber s 2006; Ravmond ct al. 2008). 

Dynamical classification: iChamber^ (|2001l) describes 
several dimensionless parameters to compare the out- 
come of simulations of accretion in the inner solar sys- 
tem to the actual planets. We adopt three, the radial 
mass concentration RMC, the angular momentum deficit 
(AMD), and the orbital spacing statistic (OSS), to clas- 
sify and compare our results. The RMC measures how 
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Fig. 8. — Plot of mass vs. semimajor axis scaled by the initial 
mass and position of the ice line, respectively, for the 6 primary 
simulation groups. The innermost planets at 0.25 to 0.60 ajj.g 
have a clear separation from the other planets. Simulations show 
a statistical mass segregation effect out to ~ 2.5 cijj^i,^^^^^. 



mass is distributed in the system and is given by 



RMC = max 



Smj.[logio(a/aj)] 



(6) 



where rrij and a,- are the masses and semimajor axes of 



Fig. 9. — Distribution of orbital eccentricities of simulated plan- 
ets for 6 = 8, 5, and b = 2 and known exoplanets (exoplanet.eu). 
Only a few of our 6 = 8,5 planets have eccentricities higher than 0.2 
whereas the observed gas giant planets span the full range of eccen- 
tricities from to nearly 1. The highly packed (6 = 2) simulations 
exhibit an eccentricity distribution much closer to observations, 
which are mostly gas giants. 



the planets in a system. A more tightly concentrated 
system will have a higher RMC. The AMD is a measure 
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TABLE 3 

Simulation Outcome Classifications 



Name 


Mtot 
(Me) 


"/ 


OSS 


cross 


RMC 


O-RMC 


AMD 


CAMD 


ST2 


10.98 


1.8 


10.32 


1.14 


13.76 


4.75 


0.0051 


0.0019 


ST3 


12.76 


2.4 


8.09 


1.81 


14.10 


5.73 


0.0095 


0.0070 


ST4 


12.12 


2.8 


7.69 


1.79 


10.79 


4.05 


0.0093 


0.0036 


LM2 


4.29 


1.8 


9.60 


1.77 


21.44 


3.87 


0.0121 


0.0076 


LM3 


4.93 


2.4 


8.84 


1.39 


17.71 


5.80 


0.0117 


0.0058 


LM4 


5.08 


2.8 


7.80 


1.96 


18.32 


3.48 


0.0152 


0.0047 


CI2 


17.16 


2.0 


8.17 


1.40 


15.61 


4.97 


0.0019 


0.0006 


CIS 


16.99 


2.5 


7.32 


1.48 


14.78 


3.48 


0.0051 


0.0021 


CI4 


16.39 


2.4 


8.00 


1.50 


13.03 


4.15 


0.0050 


0.0031 


MP2 


12.88 


1.8 


8.41 


0.98 


16.57 


3.64 


0.0038 


0.0018 


MPS 


14.22 


2.7 


7.43 


1.60 


10.76 


3.00 


0.0080 


0.0037 


MP4 


15.21 


2.9 


6.67 


0.99 


10.68 


3.36 


0.0076 


0.0041 


HPT 


26.71 


4.6 


3.85 


0.60 


9.68 


4.00 


0.0485 


0.1397 


HP9 


18.67 


4.0 


5.86 


2.19 


5.34 


4.12 


0.1560 


0.1747 


HP12 


14.89 


4.2 


6.11 


2.32 


4.74 


1.97 


0.1413 


0.1222 


LD2 


13.45 


1.9 


10.12 


0.70 


11.58 


1.58 


0.0037 


0.0013 


LD3 


13.17 


2.6 


8.27 


1.83 


10.47 


2.67 


0.0093 


0.0049 


LD4 


13.55 


3.0 


7.36 


1.99 


10.34 


6.03 


0.0093 


0.0042 



12.01 
12.05 
11.23 
11.74 
10.55 



2.4 
2.4 
3.7 
30.9 
2.7 



8.66 
7.84 
5.97 
1.11 
7.43 



2.25 
1.73 
0.83 
0.11 
0.75 



14.12 
14.78 
10.28 
8.38 
21.12 



6.01 
5.34 
1.33 
4.00 
6.92 



0.0030 
0.0031 
0.0167 
0.0521 
0.0082 



0.0050 
0.0018 
0.0096 
0.0124 
0.0038 



artifact of the detection bias for close-in planets. Our 
systems have intermediate values of RMC (8-20) that 
are poorly represented by the current catalog of known 
multi-planet systems. 



Known Wultiplanet Systems 
ST 



onfiguration at 100 Myr. 
onfiguration at 200 Myr 




of orbital excitation and is given by 



AMD 



Sj-mj 



Vaj[l - .^(1 - ej) 



cos I i 



^jlTij, 



(7) 



where ij and Cj are the inclinations and eccentricities of 
the planets in a system. The AMD measures the differ- 
ence between the angular momentum (in the z-direction) 
of a system and that of a system of identical bodies on cir- 
cular, non-inclined orbits with the same semimajor axes. 
The OSS is a measure of the mean spacing of planets and 
is given by 



OSS = 



1 



iV - 1 V a 



+ ati 



3M, 
2m 



1/4 



(8) 



where N is the number of bodies, flmax is the maximum 
semimajor axis, amin is the minimum semimajor axis, M* 
is the mass of the central star, and fh is the mean mass 
of the oligarchs. Unlike the RMC, the OSS ignores the 
mass and location of individual oligarchs and depends on 
the distance between the bodies. 

These statistics have no meaning for single-planet sys- 
tems and those cases are excluded from the calculations. 
Average values and standard deviations of RMC, AMD, 
and OSS for each set of simulations are listed in Table [3] 
along with values for a number of exoplanetary systems 
and the solar system. These data are also plotted in Fig- 
ure [TUl Although there is considerable spread in these 
statistics between the simulations, with the exception of 
the HP (b = 2) run, they occupy a region not spanned 
by known planetary systems, most of which contain gas 
giants. Our predicted systems all have OSS > 6, in con- 
trast to known expoplanet systems, but this is likely an 



Fig. 10.— Average OSS, RMC, and AMD for the 6 primary 
simulation groups alongside 12 known planetary systems with > 3 
planets as well as the solar system values. Most of the observed 
exoplanet systems shown contain at least 1 gas giant. Since incli- 
nations are measured from an invariable plane, which is not known 
or poorly defined for exoplanetary systems, we assume zero incli- 
nations for these calculations. Although the known planets cover 
a wide range of values, they are clearly very different from our 
simulations with the exception of the HP {b = 2) runs. 

Mean motion resonances: Because of migration, two 
planets may enter an MMR where the orbital periods 
are integer ratios. The resonant angle 0: 



pLi — qL2 — muji — nf^i — ru)2 — s^2 



(9) 



must librate between two values, where p, q, m, n, r, s 
are integers, L is the mean longitude, uj is the argument 
of periapsis, and n is the longitude of the ascending 
node. The subscripts 1, 2 refer to the inner and outer 
planet respectively. Outside MMR, the res onant angle 
will be unbounded (i.e., will circulate) (jElliot et al.l 
120051) . We searched the final billion years of the 180 
primary simulation sets for 1:1, 2:1, 3:2, 3:1, and 4:1 
MMRs. We required any MMR to last at least 100 kyr. 
None of our systems appear to have had a MMR in 
that interval. Even a system that appears to be in 1:1 
commensurability (see far bottom right of Figure [7|) , did 
not have a bounded resonant angle. 

Inner system: Two sets of runs (Figures II H and [T2]) in 
which mass (planets or oligarchs) was placed inside the 
ice line show that this has little effect on the dynamical 
evolution and final configuration of the outer planets. 
The innermost ice line planets were 0.1-0.2 AU further 
out at 100 Myr, and in the 103 set they had accreted an 
average of 0.3 Mq, more mass (from the inner system). 
Other effects on the outer planets were non-systematic 
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or negligible. However, the effect of the outer planets on 
the inner system were significant. Figure [TT] shows the 
configuration of the EV3 system after 100 Myr. In all 10 
runs the Earth and Venus analogs collided after 15-70 
Myr (and accreted some small bodies), forming a single 
2-3 Mq planet at a median semimajor axis of 0.81 AU 
(Figure [TT]) . In each simulation, this coincides with the 
time when the innermost of the outer planets migrates 
inside of 2.5 AU. In companion runs with Earth and 
Venus analogs but no outer planets, no such collision 
ever occurs. After the small bodies have been artificially 
removed, the configurations remain stable for at least 
1 Gyr. Figure [TH shows the 100 Myr configuration of 
systems that started with a disk of inner oligarchs rather 
than two planets. On average, 1.5 M^ of the initial 
2.2 M(^ inner disk mass has been scattered outward or 
accreted by the outer oligarchs. Amongst the 10 runs 
the largest surviving body in the inner system has a 
mass of 0.53 Mm. 
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Fig. 11.— Configuration of the EV3 set after 100 Myr. The 
planet formed from the coUision of the Earth and Venus analogs is 
shown in gray. This collision occurs in all 10 simulations, forming 
a 2 — 3 A/0 planet. The innermost planet in the system idoes not 
migrate in as far inwards as in the runs where no inner mass is 
included, but otherwise the system is similar to the ST3 set. 



Accretion onto the central star: There are two com- 
peting explanations for the observed correlation between 
high metallicity and the presence of giant planets. One 
is that higher metallicity augments the mass of solids 
in a planet-forming disk, allowing cores to form gas gi- 
ants before the gas dissipates. The other is that rocky 
material has been accreted onto the stellar photosphere 
during planet formation {Gonz alez 2006) . Because sys- 
tems without detectable gas giants are not statistically 
more metal-rich than solar, one check of our scenario 
is the amount of mass that falls onto the central star. 
An amount sufficient to significantly increase the metal- 
licity of the photosphere would conflict with observa- 



FlG. 12. — Configuration of the set starting with 130 oligarchs 
in the inner system after 100 Myr. Apoapsis and periapsis lines 
like those in Figures [7l and [TT1 are suppressed. The surviving inner 
system oligarchs are shown in gray. In most cases the inner system 
oligarchs were thrown onto high inclination orbits outside of the 
inner system or accreted by an inward migrating oligarch. There 
is no obvious pattern to the distribution of the (initially) inner 
system oligarchs. 

tions. We assume that a solar ma ss star had a convec- 
tive region of 0.02 M^o at 1 Gyr, (jAsplund et al.l 120091 : 
iD'Antona fc Mazzitellil 119941 ) and that solids have the 
compos ition of carbonaceous cho ndrites (Fe is 20% by 
mass) (jGrevesse fc Anderill989| ). In our 180 primary 
simulations without mass in the inner system, the aver- 
age mass accreted onto the parent star is 0.9 ± 0.8 M®. 
For the EV3 and 103 sets, the average is 0.3 ± 0.1 M® 
and 3 ± 1 Mq respectively. The highest value in any run 
is 4.8 M0, corresponding to ~ 1 M^ of iron. The cor- 
responding increase in [Fe/H] is no more than 0.06 dex 
and more typically ~ 0.()2 dex. The actual metallicity en- 
hancement is likely to be smaller because the convective 
zone of solar-mass stars is much larg er at t < 30 My r 
when much of the mass is accreted (jFord et al.l 119991 ). 
Thus our scenario does not conflict with observations. 



3.3. Sensitivity Runs 

The TP3 runs contain twice as many small bodies 
(1000) than the other runs, but produce systems with 
the same mean number of planets, and similar values of 
Mtot, OSS, RMC, AMD compared to the ST3 set (Table 
[21). The only major difference was that there was less 
variation between the 5 Gyr configurations produced by 
the TP3 runs. The standard deviation of AMD, OSS, 
and RMC are all lower, presumably because random fiuc- 
tuations are reduced with a larger number of small bod- 
ies. We conclude that in most of our simulations, the 
between-run variability is exaggerated due to the use of 
a finite number of small bodies. 

Simulations in which the small bodies were removed at 
2 Gyr as opposed to 1 Gyr did not result in significantly 
different systems. Most systems, including both oligarchs 
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and small bodies, achieve a degree of stability earlier than 
1 Gyr. Runs where the initial masses, position, and Hill 
spacing of the oligarchs vary slightly (Table [Ij did not 
produce significantly different systems. 

Simulations run with higher eccentricities and inclina- 
tions for small bodies nearby the oligarchs showed slower 
initial mass growth than the standard set. At 100 Myr, 
the HE3 set had ~ 4 M^ less mass in oligarchs than the 
ST3 set. By 250 Myr, the difference was only ~ 2 M®. 
This result supports our expectation that because the 
Safronov numbers are low, oligarchs ultimately accrete 
all planetesimals in their zone. Inner planet migration is 
minimally affected by the higher eccentricities and incli- 
nations. At 100 Myr the innermost planet is, on average, 
< 0.2 AU further out in the IIE3 set than in the ST3 set. 
At 250 Myr, the difference is negligible. 

4. PROSPECTS FOR DETECTION 

Our simulations are useful to the extent they can make 
testable predictions. We investigate the prospect of de- 
tecting the 534 planets predicted by our 180 primary 
simulation sets. We divide the expected detections up 
by initial conditions to determine which initial condi- 
tions were most accurate when detections are made. 
Figure [13] shows their masses and semimajor axes rela- 
tive to the detection domains of Doppler radial velocity, 
ground-based microlensing techniques, the NASA Ke- 
pler mission, and the proposed SIM-Lite. Until there 
are substantial improvements in sensitivity and stabil- 
ity (JEggenberger & Udry 2010), Doppler is unlikely to 
detect any of the predicted planets. The other three 
techniques will be able to detect at least some of these 
objects. 

4.1. Microlensing 

Microlensing is currently the only ground-based de- 
tection method that is sensitive to the planets pre- 
dicted by our simulations. Microlensing is most sen- 
sitive to planets with projected separations near the 
Einstein radii i?E of their primaries, corresponding to 
Re, ^ 3.5 AU(M,/M0)^/2 for typical lens and source 
distances. Thus for a typical primary mass in current 
surveys of ~ 0.5 Mq (jGould et al.ll20ld[ ). the sensitivity 
of microlensing peaks for planets with semimajor axes 
^ 3 AU. Current microlensing surveys can detect plan- 
ets with mass > 3 M@ with separations within a fac- 
tor of a few of this distance. Indeed, several of the mi- 
crolensing planets detected to date have masses in the 
range 3 — 15 Mm and projected sep arations of 1 — 3 AU 
(iBeaulieu et all 120061 fGould et all l2006t iBennett et all 
120081 ). and thus may be analogs to our simulated sys- 
tems. Interestingly, for the planet OGLE-2005-BLG- 
169Lb with mass ^1 3 Mq and projected separation ^ 
2.7 AU, [Gould et aj ] (|200al exclude additional Jupiter- 
mass planets within the range of projected separations of 
0.5 — 15 AU; indicating that this may indeed by a system 
without gas giants. 

We can estimate the expected number of microlensing 
detections one would expect, assuming that 60% of stars 
have systems such as those we simulate. Using the stan- 
dard set, for each system we randomly choose a primary 
lens mass according to an event rate distribution 
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Fig. 13. — Final mass and semimajor axis of all planets from 
the primary 6 sets and the detection domains of 4 planet-finding 
techniques. The Doppler range is set by a Doppler amplitude of 
K = 3 ms~^. A Kepler detection requires observation of at least 2 
(rather than the usual 3) transits, and we assume a mission lifetime 
of 3.5 years, so planets with period P > 1.75 yr will not be detected. 
Microlensing is most sensitive in the 1.5-6 AU range, where the 
planet detection probability is at least 1% per microlensing event. 
The SIM- Lite range is set by a detection probability > 85% (see 
Equation l|17p in the text). Although SIM might be able to detect 
planets with P greater than the lifetime of the mission (~ 5 yr) we 
conservatively exclude such planets. 

for a mass function dN/ dlogM (x Af""^^. We adopt 
a = 0.2 and restrict our primary masses to the range 
0.05— 1 Mq, with an average primary mass of ~ 0.5 M©. 
We assume the planets in the system are coplanar and 
draw a random inclination i for the system distributed 
as cosi. Then, for each planet, we compute its mass 
ratio and projected separation, drawing a random or- 
bital phase for each planet, ignoring the (small) ef- 
fects of non-zero eccentricities. We then scale the pro- 
jected separation to the Einstein radius, assuming i?E — 
3.5 A\J{M^./MqY^^. Finally, we determine which of the 
plane ts in the system are detected in each of the 13 events 
in the lGould et al.l (|2010t ) sample, noting instances when 
multiple planets are detected. We repeat this for all of 
the simulated syst ems and for 5000 Monte Carlo tri- 
als. We find that iGould et all (|2010[ ) should have de- 
tected ~ 1.7 planets, with an expected mean mass ratio 
of ^ 10""*^, and maxim um mass ratio of ~ 10~'^'^. In 
fact. IGould et a l. (2010") found one system with mass ra- 
tio ~ lO--*-! (OGLE-2005-BLG-169Lb), and two systems 
with mass ratio lO"'^'^, consistent with our scenario. 

Figure [HI shows the observed curn ulative distributions 
of mass ratios from the lGould et al.l ((lOlO) sample, com- 
pared to the expected distributions for a scenario in 
which 60% of stars host planets with the properties of 
the standard simulation set, and 30% host four giant 
planets with the masses and semimajor axes of the solar 
system. The remaining 10% of giant planet systems host 
close-in planets currently undetectable by microlensing. 
The number of expected detections and the distribution 
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of mass ratios are both broadly consistent with the ob- 
served sample of events. We conclude that this scenario 
is consistent with all available constraints. 
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Fig. 14. — The solid line shows the cumulative distribution of 
mass ratios for the six detected planets in the four year sample 
of 13 microl ensing events monitored by the /iFUN collaboration 
H Gould et al . 2010). The dotted line shows the cumulative distri- 
bution of mass ratios predicted for this sample, based on the detec- 
tion efficiencies of the monitored events, and assuming a model in 
which 30% of stars have four giant planets with masses and semi- 
major axes equal to Jupiter, Saturn, Uranus, and Neptune, and 
60% of stars have systems of planets predicted by our standard 
simulation set. These predictions assume a power-law distribution 
of primary masses, with a mean mass of ~ 0.5Mq. 



What are the prospects for detecting analogs to the 
systems we have simulated in future microlensing sur- 
veys? We consider the two classes of microlensing sur- 
veys that are likely to take place over the next 10 years: 
alert and follow-up monitoring of high-magnification 
events similar to that already being conducted, and 
"next-generation" surveys in which thousands of low- 
magnification events are detected and simultaneously 
monitored with the ^10 minute cadence needed to de- 
tect Earth-mass planets using an array of 1 — 2m tele- 
scope s with wide field-of-view cameras. See (jGaudi et al.l 
I2009D for further discussion of these two channels. 

For the high-magnification event channel, we follow 
the method outlined above to simulate the number of 
expected detections, except we assume that 20 events 
per year with maximum magnification > 100 are densely 
monitored during each peak. This r epresents a f a ctor o f 
~ 6 improvement over the rate in iGould et al.l ()2010[ ). 
which should be realizable with the expected better pre- 
diction of high-magnification events, increased number 
of alerts, and decrease in the maximum m agnification 
threshold from 200 to 100 (|Gould et al.|[2010l ). We adopt 
the analytic detec tion sensitivity estimate discussed in 
(jGould et al.|[20l0l ). assuming r? = 0.32 and £, = 100. The 
results are shown in Table |4l for the six primary simula- 
tion sets. We expect an average of 4.6 planet detections 



per year (for the standard set), with roughly one detec- 
tion of a multiple-planet system per year. For the HP 
simulation set, we find that there is a significant chance 
(0.17 per year) of detecting as many as four planets in 
the same event, whereas these probabilities are generally 
substantially smaller (<0.03 per year) for the other sim- 
ulations. This indicates it may be possible to distinguish 
between the various input assumptions of the simulations 
using observations of multiple planet systems. 

For the low-magnification events detected in next- 
generation surveys, we use the unpublished simulation 
code of Gaudi, Han, and Gould. This code simulates en- 
sembles of planetary microlensing events and estimates 
detection rates for a given input value of the mass and 
semimajor axis of the planet. The simulated light curves 
account for the effects of weather, variable seeing, moon 
and sky background, and the finite size of the source 
star. We assume parameters similar to that expected 
for the funded Korean Microlensing Telescope Network 
next-generation microlensing survey: three l.Gm tele- 
scopes with 4 deg^ cameras located in Australia, Chile, 
and South Africa (C. Han, pers. communication). The 
host lenses are drawn from a model of the Galactic pop- 
ulation of lenses and sources that matches available con- 
straints pan fc Gouldlll995l[200l . The resulting detec- 
tion probability for each of the planets in the 180 simu- 
lations is shown in Figure I15[ and the expected number 
of detections per year are listed in Table HI We predict 
that next-generation surveys should detect ^ 22 planets 
in low-magnification events per year (for the standard 
set), assuming that 60% of all stars host planetary sys- 
tems such as those we simulate. These detections are 
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Fig. 15. — Probability of detecting the planets from our 180 
primary simulation sets in a next-generation, ground-based mi- 
crolensing survey, consisting of three 1.6m telescopes with large 
FOV cameras located in Chile, South Africa, and Australia. The 
higher detection probabilities near ~ 2 — 3AU are caused by their 
proximity to the Einstein ring and t he t endency for closer in plan- 
ets to have higher masses (Section 13.211 . Assuming 60% of stars 
systems analogous to those in our primary simulation sets, such 
a next-generation ground-based microlensing survey would detect 
~ 22 planets per year (see Table [4]|. 
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TABLE 4 
Detection Statistics 



Name 


Kepler^ 


SIM^ 


Microlensin^ 
(low-mag) (high-ma 
Total M-P 


Sys^ 


TTV^ 
(min) 


Transit Prob^ 


ST 


128.9 


39.8 


22.4 


4.63 




1.07 


34.3 


0.30 


LM 


7.0 


34.3 


13.2 


2.15 




0.52 


73.2 


0.19 


CI 


344.9 


55.3 


23.7 


5.86 




1.66 


14.8 


0.47 


MP 


105.7 


39.4 


24.9 


5.58 




1.48 


32.6 


0.29 


HP 


42.8 


32.3 


22.8 


5.98 




1.26 


54.2 


0.21 


LD 


119.2 


38.6 


22.4 


4.89 




1.13 


38.1 


0.28 


Average 


124.8 


39.9 


21.1 


4.62 




1.08 


41.2 


0.28 



[jHere, a Kepler detection counts if > 2 transits are observed over the 3.5 yr Kepler mission. We assume ^ 60% of stars have 
systems similar to those in a given simulation set. 

[^Number of detections by SIM-Lite assuming 64 target stars and that ~ 60% of stars have systems similar to those of a given 
simulation set. 

jThe number of microlensing detections per year assuming ~ 60% of stars have systems similar to those of a given simulation set. 
jMedian transit timing variation for the innermost planets in a given simulation set. 
jMedian transit probability for innermost planet in a given simulation set. 
jNumber of systems with more than one planet detected in a single microlensing event. 



in addition to those found in high-magnification events. 
Muhiple-planet systems wiU be rare (detection probabil- 
ities of < 0.1%) in these low-magnification events. 

While ground-based surveys are relatively insensitive 
to the low-mass, large semimajor axis planets we typi- 
cally find in our s imulated systems, a space-based mi- 
crole n sing survey (iBennett fc Rhiell2002t iBennett et alJ 
120091 : iBeaulieu et alJ l2010[ ) would be exquisitely sensi- 
tive to these bodies (and essentially all of the planets 
we find in our simulations). In particular, a space-based 
microlensing survey would detect the most distant plan- 
ets with a > 15 AU as isolated, short time scale event s 
without the signature of the host star (jHan et alJl2005j ). 



Ra 



P 



3.5 yrs 



1/6 



-^Q-0.2(y-12) 



(11) 



where i?* and M* are the radius and mass of the star, 
Tp is th e radius of the planet, and P is the period of the 
planet (|Yee fc Gaudi|[200l . The radius of a 10 Afgj body 
composed of equal p arts water ice and ro ck/metal is pre- 
dicted to be 2.3 i?0 (jGrasset et alJl2009( ). At the median 
semimajor axis of our innermost planets (a = 1.66 AU, 
P = 2.14 yr) a planet orbiting a F = 12, Sun-like star 
would produce a transit with a depth of 0.4 mmag and 
a S/N of - 50. The depth and S/N will be larger if the 
planet has a thick atmosphere. 



4.2. Kepler 

The Kepler spacecraft was successfully launched on 
2009 March 6 and is continuously monitoring ~ 10^ F- 
to K-type stars with the primary objective of discov- 
ering transiti ng Earth-mass pl anets on 1 AU (1 yr pe- 
riod) orbits (jKoch et al.ll2010l ). although the detection 
of many planets on s horter-period orbits is expected, 
e.g., iSelsis et all ()2007| ). Three transits will be required 
to confirm a plane t; hence the 3 . 5 yr nominal mission 
lifetime. However, lYee fc Gaudil (|2008| ) point out that 
Kepler should detect one or two transits by planets on 
more distant orbits. The innermost planets in our 180 
primary simulation runs have a median a = 1.66 AU 
(P ~ 2.14 yr), making it possible that two (but usually 
not three) transits would be observed, geometry permit- 
ting. We calculated the expected number of such planets 
that Kepler will detect transiting at least twice using 
Equations (2) and (4) from Yee & Gaudi (2008), assum- 
ing that 60% of all solar-type stars have such systems, 
and ignoring the effect of eccentricity. (The median ec- 
centricity is . 08). We use Kepler's precision given in 
IJenkins et al.l (I2010D and the ch a racter istics of Kepler's 
target stars from iBatalha et al.l ()2010l ) . The predicted 
number of transit detections is 129 for the standard set, 
and will be larger if ice lines are located closer to stars 
(See Table H]). Around a solar mass and radius star Ke- 
pler's detects a transit with a signal-to-noise ratio (S/N) 
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Fig. 16. — Transit timing variations vs. transit probability of the 
innermost planets predicted by our simulations. The TTV is the 
absolute difference in the time of transit center between successive 
transits, averaged over successive transits. The estimated TTV 
precision of Kepler (Ic) is 4 minutes in short cadence mode. TTV 
scales as orbital period P, and transit probability as P~ ' , so that 
TTV ~ (transit probability) "^Z^. 
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In addition to the orbital period P (the interval be- 
tween transits) it may be possible to limit the orbita l 
eccentricity of such transiting planets (jFord et al.ll2008[ ). 
To the extent that limb darkening effects can be ac- 
counted for, the orbital eccentricity can be constrained 
independently of the transit impact parameter by mea- 
suring the duration of both the transit T and the ingress 
and egress phases r. The following relation holds: 



(1 4- esinujy 



Tt 



3P 



2/3 



(12) 



where p* is the density of the star, and a; is the ar- 
gument of the periapsis (jSeager fc Mallen-Ornela£|[2003l : 
iTinglev fc Sacket3l2005D . To the first order in eccentric- 
ity, the term on the left is 1 — 2e sin a; and therefore only a 
minimum eccentricity can be established independently 
of the longitude of periapsis. However, the distribution 
of w in a population of transiting systems must be uni- 
form over — 27r and thus the statistical distribution of e 
can be inferred, p* is well established from astrophysical 
theory but will not be measured directly except in those 
cases where there is a second transiting planet either on 
a shorter period orbit that is accessible to the Doppler 
technique (jSozzetti et a l."2007') or via thei r mutual tran- 
sit timing perturbations (Holnia n fc Murr av 2005). 

The presence of an additional, outer planet of 
mass M2 on an orbit with semimajor axis 02 can 
be inferred th rough variations in the interval be- 
tween transits ("Mirald a-Escudel [200l lAgol et al.l 120051: 
iHolman fc Mur rav 2005]). The standard deviation of the 
transit time due to perturbations from the outer planet 
is approximately 



STr, 



Pi 02 M2 

23/27r ai M» 



(13) 



for nearly circular orbits (jAgol et al.ll2005[ ). This will be 
typically 3 x 10~^Pi, or ~ 5 min. If the second planet 
is near a mean-motion commensurability, the variation 
can be 1-3 orders of magnitude larger. We calculated 
the STc of the inner planet induced by the other planets 
in each of the 180 systems. These calculations used the 
orbited elements computed by Mercury for 100 successive 
orbits at 5 Gyr. The linear perturbations of the orbital 
equations are: 



where 



and 



6Tc^5to + ^S^i+^SP, 

Ztt Zn 



5fi = (1 — e cos ri)6ri — (sin ?/) 6e, 



(14) 
(15) 



6rj = — 



cos uj (l — e cos rj) Suj 
( 1 + e) sin T] — e sin rj cos 77 

[1 — (1 -I- e) cos 77] Se 
(1 -I- e) sinry — esinijcosr]' 



(16) 



where to is the time of periapsis passage, fi is the mean 
anomaly, and rj is the eccentric anomaly. Figure [12] 
plots the 6Tc versus transit probability for the inner- 
most planets in the 180 simulation sets. Table |4] has 



the 6Tc and transit probability values for each simula- 
tion set. Detection of TTV with Kepler obviously re- 
quires a third transit, and thus an extended mission, as 
well as read-out in short-cadence (59 s) mode to capture 
the ingress or egress (about 20 min in duration). The 
S/N of the transit detection over the ingress and egress 
is about 5, and the la precision of the timing (using both 
ingress and egress) is about 3 min. Thus additional ob- 
servations of transits by the innermost planet by Kepler 
should be sufficient to reveal the presence of outer plan- 
ets like those predicted by our simulations. Observations 
from the ground have achieved ^ 0.5 mm ag precision 
pohnson et al.ll2009HSo"uthworth et 'alll2009[ ). raising the 
possibility that ground-based follow-up might also reveal 
such variation. 

4.3. Space Interferometry Mission (SIM-Lite) 

SIM-Lite is an astrometric interferometer mission that 
will achieve sub-microarcsecond precision per "visit" 
and should be capable of detecting Earth-mass plan- 
ets in the habitable zone of nearby [d < 30 pc) stars 
(jShao fc Ne mati 2009). During a nominal program that 
consumes 40% of a 5 yr mission lifetime, 64 target stars 
can each be visited about 200 times. The total S/N in 
N visits is 

S/N-F.^^^, (17) 

where rup is the mass of the planet, a is the measurement 
error in arc-seconds, D is the distance in pc, and the 
dimensionless factor F is 



(l + cos2i) l-e2 3- ^-1 cos 2a; 



1 + cos2 i 



(18) 
(K. Mogren fc B. S. Gaudi, i preparation, J. Catanzarite, 
private communication). Although planets with orbital 
periods longer than the mission lifetime might be de- 
tectable, we conservatively assume that this is not the 
case. Equation P^ is averaged over an isotropic distri- 
bution of values for i and w, and the S/N is calculated 
for the 64 target stars for SIM-lite provided by J. Catan- 
zarite. We ado pt a detection criteria o f S/N > 5.8 for 
a single planet (jCatanzarite et all I2006D . We find that 
96% of the innermost planets will be detected. If we 
assume ~ 60% of stars have systems like those in our 
simulations and ignore planets with periods greater than 
the lifetime of the mission, we predict that SIM-lite will 
find 44 planets like those in our simulations (see Table 
U). The presence of multiple planets may make disam- 
biguation of orbital parameters difficult, but detection 
of addit ional planet s in these systems is clearly possi- 
ble. See lFb^ (I2OO6I) and lGould((l200l for more robust 
analyses of this problem. 

5. DISCUSSION 

5.1. Summary 

Our simulations predict the evolution, final configura- 
tion and detection of systems of icy planets lacking gas 
giants. If, as observations indicate, all solar-mass stars 
are born with disks, but only a minority (^40%) form 
giant planets, our predictions may describe the hitherto 
"invisible" majority of outer planetary systems. 
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Because planet formation is so poorly understood, and 
there are few constraints on our simulations, we explore 
a wide range of initial conditions. We assume the canon- 
ical theory of planet formation in which a small number 
of oligarchs grow from a disk of much smaller planetes- 
imals. Motivated by the core accretion theory of giant 
planet formation, we also assume that oligarchic growth 
was accelerated in a region of enhanced surface density 
of solids beyond the ice line, a hypothetical point in the 
planet-forming disk where water ice condenses. We do 
not simulate the formation of these oligarchs, only their 
subsequent growth and dynamics after the gas has disap- 
peared from the disk. We assume that the longer orbital 
periods and lower surface density has stymied oligarch 
formation further out in the disk. The residual disk is 
represented by a finite number (typically 500) of small 
bodies. For increased computational efficiency, we place 
mass interior to the ice line in only a few realizations to 
investigate its effect. 

We vary the total mass in the ice line, the location 
of the ice line, the initial number of oligarchs, and their 
spacing. We evolve each system for 5 Gyr; after the first 
Gyr the small bodies are removed. These systems are 
highly chaotic in the first 10 Myr, but usually become 
stable well before 1 Gyr, and removal of small bodies has 
no significant effect on the oligarchs' subsequent evolu- 
tion. We describe the evolution of these systems in terms 
of the time to clear the oligarch zone of small bodies, 
the number of MMR crossings, the number of oligarch 
"swaps" , the migration of the innermost oligarch, and 
the frequency of oligarch ejection (Tabled. We describe 
the final configuration of the systems at 5 Gyr using three 
statistics: the OSS, RMC, and AMD (Table[3l). In a lim- 
ited set of observations, we place mass interior to the ice 
line, either pairs of Earth Venus analogs, or smaller oli- 
garchs, and we investigate its effect on the evolution of 
the outer system, as well as its own fate. 

5.2. Major conclusions and implications 

In the vast majority (169/180) of our primary runs, 
and across all initial conditions we investigate, an oli- 
garch migrates interior to the ice line, settling to be- 
tween 25% and 60% of the ice line distance in about 10 
Myr and growing into a planet with a median mass of 
0.23 Mice. We call this object the innermost migrated 
planet, or IMP. In 123 of the 169 primary runs with an 
IMP, the IMP is the most massive planet. IMPs are 
clearly distinguishable from the other planets (Figures [7] 
and [S]) . The migration is a result of exchange of angu- 
lar momentum between the IMP and the other, exterior 
oligarchs: 5 of the 11 runs that did not produce an IMP 
contain only a single planet at 5 Gyr. The migration 
is significant because, unlike with gas giants, the plan- 
ets in our simulations have low masses compared to the 
residual disk mass. The existence of mass in the inner 
system only slightly affects the final position and mass 
of IMP, but the converse is not true (see below). IMPs 
may be the visible representatives of an otherwise "invis- 
ible" majority: The common occurrence, relatively high 
mass, and small semimajor axis of IMPs make them emi- 
nently detectable by microlensing, transits (with Kepler), 
and astrometry (with SIM-Lite), but not yet by current 
Doppler capabilities. 

Ground-based microlensing is currently capable of de- 



tecting planets as small as '-^ 3 Af^ at separations of 
1.5-3 AU, and indeed several planets with (uncertain) 
masses between a few M^ and one or two Neptunes have 
been found in this distance range (iBeaulieu et "all 120061 : 
iGould et al .' 2006; 'Benn ett et alJl2008l : iSumi et al.ll2010( ). 
These few detections may represent only the tip of the 
IMP-berg: all available constraints on the frequency of 
gas-giant and lower-mass planets from current radial ve- 
locity and microlensing surveys are consistent with the 
scenario that the minority of stars host gas giants, and 
that at least ~ 60 % of stars host sys t ems such as thos e 
we have simulated (IGould et al.ll2006l : iSumi et al.ll2010D . 
Future microlensing surveys will provide a definitive sta- 
tistical measurement or upper limit on the frequency of 
systems like those predicted here. If 60% of stars indeed 
host systems similar to those we simulate, we estimate 
that ne xt-generation ground-based m icrolensi ng surveys 
(IGould et al. 2007: Gould 2008: Gaudi et all [20091 wiU 
detect ~ 26 planets per year, including a handful of 
multiple-planet systems. A space-based microlensing 
survey would be se nsitive to essentially ah of the planets 
wc have simulated ("Ben nett fc Rhiell2002l : iBennett et al.l 
2009; Bcaulicu et al. 201^^ 

Assuming an ice-rock composition, all IMPs predicted 
here would produce a transit sufficiently deep to be de- 
tected by Kepler. 83% have periods less than the space- 
craft's 3.5 yr mission. If IMPs are present around 60% of 
solar- type stars, we predict that Kepler will detect ~129 
of them with two or more transits. Observations of addi- 
tional transits in high cadence (1 min resolution) mode in 
an extended Kepler mission could reveal additional, ex- 
terior planets through the variation of the timing of tran- 
sits. Direct calculations show variations of 20-90 min in 
our predicted systems. Finally, SIM-Lite should be able 
detect 96% of IMPs. 

The planets we predict have, statistically, very differ- 
ent orbital, mass, and eccentricity distributions than the 
giant planets discovered to date by Doppler surveys. The 
orbital eccentricities of our predicted planets are signif- 
icantly lower than in known exoplanetary systems. Our 
results agree qualitatively with the simulations of sys - 
tems of Neptune-size planets by iRavmond et al.l (|2010f ) . 
The exceptions are, unsurprisingly, our sets of simula- 
tions with dynamically overpacked oligarchs, whose or- 
bital eccentricities at 5 Gyr resemble the observed exo- 
planet distribution. Kepler observations of the duration 
of transit ingress and/or egress will offer limited con- 
straints on the distribution of eccentricities of IMPs, if 
they are sufficiently numerous. 

If many more IMPs are found than predicted here, and 
they are closer to their parent star, then a possible ex- 
planation suggested by our simulations is that the ice 
line is often much closer than 5 AU from the parent star. 
Conversely, if the combination of microlensing, Kepler, 
and SIM-Lite fail to discover a population of IMPs, one 
or more assumptions in our scenario is false. The most 
likely suspect would be the assumption of a significant 
concentration of solids at or immediately beyond the ice 
line. This would have ramifications for the core accretion 
theory of giant planet formation. 

Intriguingly, inner systems of two dominant planets 
are not stable in our scenario. In all simulations with 
Earth and Venus analogs, the two bodies collided, form- 
ing a single body at ~ 0.8 AU. No contradiction with 
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the solar system is engendered because it contains giant 
planets. Such a conglomerate would induce a bary cen- 
ter motion of 0.2 m s~^ (260 d period) which may be 
detect able by future ultra- high precision Doppler moni- 
toring (jPepe fc Lovisl[2008l ) . On the other hand, if planet 
formation in the inner system has progressed only to the 
giant impact (oligarch-dominated) phase by the time the 
IMP migrates inward, the IMP will clear most of this 
mass, leaving only small (< 0.3 M^) bodies. 

If disruption of the inner system does occur, the IMP 
is left as the only detectable planet near, but exterior 
to, the nor ninal habitable zone of an Earth "twin" (0.95- 
1.37 AU) (jKasting et al.l [19931 ). However, it is expected 
that surface temperatures will be higher on more mas- 
sive planets with thicker atmospheres, such as could be 
the case for the IMP. Given that IMPs have large quanti- 
ties of water, it then follows that IMP-like planets could 
be the most numerous type of habitable planet in the 
universe. 

5.3. Limitations of our simulations 

We have ignored Type I migration in our simulation. 
However, the magnitude and even the sig n of Type I mi- 
gration is not yet clear (Li ct al. 2009; Muto k. Inutsukal 
I2OO9I: [()gih"ara fc Idal 120(11 lYu et al.ll2010l ). Our results 
are best seen in the context of being an end-member of a 
larger suite of scenarios in whi ch Type I migratio n plays 
a role to a varying degree, c.f. Ilda fc Liiil ()2008af) . 

Planetesimals will fragment (rather than accrete) if the 
collisional energy exceeds the strength of the colliders, 
and the production of smaller fragments can eventually 
produce a collisional "cascade" whose ultimate product 
is micron-sized dust which will be swept from t he disk 
by ra diation forces or coupling to the gas disk (jWvatti 
l2008t ). After the gas disk disappears, the orbital ec- 
centricities and inclinations of planetesimals are no long 
damped, and they will be excited by the oligarchs, which 
are growing by accretion of planetesimals. Thus fragmen- 
tation will compete with oligarch accretion and may limit 
oligarch growth to some maximum mass. In a collisional 
cascade, most of the mass will be in the largest planetesi- 
mals and thus it is their lifetime that will set the balance 
between accretion and fragmentation. The strength of 
these will be set by gravity and will depend on on size as 
s'^/^ ([Krivov et al.l [20051 ) . and thus the critical collision 
speed for fragmenta tion will scale as s '^^^, e.g., 1 km s^^ 
for a 100 km body (jLohne et al.ll2008| ). The equilibrium 
velocities of the planetesimals will scale with those of 
the oligarchs by the ratio of mass surface densities in 

the respective populations ( ^^ ) where n « 0.25 — 0.5 



(e.g. JGoldreich et al.l (J2004D '). Thus as planetesimals are 
accreted and oligarchs grow, the collision speeds of the 
former will increase until the largest planetesimals suf- 
fer destructive collisions, after which that process com- 
petes effectively with accretion for mass. Unfortunately, 
neither the size of the largest planetesimals (which may 
ultimately derive from physics such as two-stream insta- 
bilities (jJohansen et al.ll2007( ). However, once E/ ~ Sp, 
Vp will become comparable to the escape speed of the 
oligarchs (~ 10 km s~^) and probably well above the 
threshold for fragmentation. Thus a crude upper limit 
on the effect of fragmentation is the mass of the oligarchs 



when the surface density of planetesimals fall below the 
surface density of oligarchs. In the standard set of sim- 
ulations, this is usually after ^ 100 Myr. Direct evalu- 
ation of the RMS encounter velocity of planetesimals in 
our standard sets show that this kinetic energy increases 
with time until it reaches 10% of the orbital kinetic en- 
ergy (i.e., 11 ~ 4 km s^^) in ^ 100 Myr. Oligarchs have 
accreted ~ 80% of their mass by this time. Orbital mi- 
gration, which depends on the mass surface density of 
the planetesimals, and not their mass distribution, will 
be affected less. 

Our scenario assumes a significant concentration of 
mass within 1 AU of the ice line, with a total amount 
sufficient to produce least one and as many as four of the 
supposed 5-10 M® cores of the outer planets in our solar 
system. If the total mass in the ice line was less, or it 
was less concentrated then we assume, the result would 
be smaller bodies, and the amount of inward migration 
by the innermost object would be less. 

Our systems are all orbiting solar-mass stars, whereas 
those s urveyed by Kepler comprise a range of stellar 
masses (jBatalha et al.l[2010l ). and the preponderance of 
microlensing stars are lower mass M dwarfs (< O.5M0). 
Observations suggest a correlation between stellar mass 
and giant planet frequen cy, at least on detectable or- 
bits ([Johnson et al.ll2007D. and in line with som e theo- 
retical expectations ([Kennedv fc Kenvonll2008b[ ). At a 
given ice line distance, the orbital time scale is longer 
and the Safronov number scales inversely with stellar 
mass. We thus expect those systems around M dwarfs 
to develop more slowly and scattering to be more effi- 
cient relative to accretion. A prediction of the latter is a 
higher mean orbital eccentricity among the planets than 
the values found here. However, a correlation between 
disk mass and stellar mass, if one exi sts (JNatta et al] 
'20001 lEisner et al.l[200llVorobyovl[2009[ ). along with our 
finding that the mass of the IMP approximately scales 
with the mass in the ice line, would partially offset this 
effect . Moreover, th e ice line itself may be closer to the 
star ([Kennedy fc K envon 2008b), and thus both the or- 
bital time scale and Safronov number will depend only 
weakly on stellar mass. 

Given these unresolved issues, our findings should be 
considered a series of predictions of one class of plan- 
ets that could be (and perhaps is being) discovered by 
microlensing, Kepler, and a future SIM-Lite mission. If 
gas giant-containing systems are indeed in the minority, 
then systems of icy Earth-to-Neptune-mass planets may 
be the hitherto undiscovered majority of planet systems, 
and their innermost members - the IMPs - could be one 
of the most common abodes for life in the universe. As- 
suming the continued success of microlensing surveys and 
planet-finding missions like Kepler, we shall soon know 
the answer. 
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